# wczytanie pakietu
library("terra")
W pierwszym kroku musimy stworzyć listę plików (rastrów), które
zamierzamy wczytać. W tym celu możemy wykorzystać funkcję
list.files(), która jako argument przyjmuje ścieżkę do
folderu z plikami. Oprócz tego musimy wskazać jaki rodzaj plików chcemy
wczytać (pattern = "\\.TIF$") oraz zwrócić pełne ścieżki do
plików (full.names = TRUE).
# listowanie plików z katalogu
files = list.files("dane/landsat", pattern = "\\.TIF$", full.names = TRUE)
files
## [1] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1.TIF"
## [2] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2.TIF"
## [3] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3.TIF"
## [4] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B4.TIF"
## [5] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B5.TIF"
## [6] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B6.TIF"
## [7] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B7.TIF"
Kiedy utworzyliśmy już listę plików, możemy je wczytać przy pomocy
funkcji rast() z pakietu terra i następnie
wyświetlić metadane.
# wczytanie danych rastrowych
landsat = rast(files)
landsat # odwołanie się do obiektu wyświetla metadane
## class : SpatRaster
## dimensions : 1853, 2846, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 602535, 687915, 5742525, 5798115 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633)
## sources : LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1.TIF
## LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2.TIF
## LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3.TIF
## ... and 4 more source(s)
## names : LC08_~SR_B1, LC08_~SR_B2, LC08_~SR_B3, LC08_~SR_B4, LC08_~SR_B5, LC08_~SR_B6, ...
## min values : 21, 32, 1321, 1924, 6827, 6982, ...
## max values : 60441, 61519, 65451, 65192, 63873, 65454, ...
Możemy również skrócić lub zmienić nazwy kanałów spektralnych. Przed tą operacją należy się upewnić czy kanały zostały wczytane w prawidłowej kolejności.
names(landsat) # nazwy oryginalne
## [1] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1"
## [2] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2"
## [3] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3"
## [4] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B4"
## [5] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B5"
## [6] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B6"
## [7] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B7"
names(landsat) = paste0("B", 1:7) # skrócenie nazw
names(landsat) # nowe nazwy
## [1] "B1" "B2" "B3" "B4" "B5" "B6" "B7"
# zamiana nazw
# names(landsat) = c("Ultra Blue", "Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2")
Wczytanie danych wektorowych odbywa się w analogiczny sposób za
pomocą funkcji vect().
# wczytanie danych wektorowych
poly = vect("dane/powiat_sremski.gpkg")
poly
## class : SpatVector
## geometry : polygons
## dimensions : 1, 1 (geometries, attributes)
## extent : 622510.6, 659710.4, 5754795, 5784772 (xmin, xmax, ymin, ymax)
## source : powiat_sremski.gpkg
## coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633)
## names : TERYT
## type : <chr>
## values : 3026
Teraz możemy przygotować prostą wizualizację na przykładzie kanału bliskiej podczerwieni (near infrared; B5) oraz poligonu.
# wizualizacja
plot(landsat[[5]]) # alternatywnie: plot(landsat[["B5"]])
plot(poly, add = TRUE)
Zasięg naszego obszaru analizy ograniczony jest do powiatu
śremskiego, natomiast scena satelitarna ma o wiele większy zasięg. W
takiej sytuacji możemy dociąć rastry, dzięki czemu ich dalsze
przetwarzanie będzie szybsze, a wynik obliczeń zajmie mniej miejsca na
dysku. Do docinania rastrów służy funkcja crop() i jako
argumenty musimy podać raster oraz wektor.
landsat = crop(landsat, poly)
plot(landsat[[5]])
plot(poly, add = TRUE)
Obszar rastra zmniejszył się. Jednak możemy zauważyć, że poza
poligonem wartości nie zostały usunięte. Wynika to z faktu, że obraz
zawsze docinany jest do obwiedni (bounding box), a wartości
poza obiektem/poligonem w rzeczywistości są maskowane (tj. są oznaczane
jako brakujące wartości). Aby zamaskować piksele poza poligonem należy
użyć funkcji mask().
landsat = mask(landsat, poly)
plot(landsat[[5]])
plot(poly, add = TRUE)
Tą operację można przeprowadzić również w jednej linii kodu używając
argumentu mask = TRUE w funkcji crop().
crop(landsat, poly, mask = TRUE)
W następnym kroku możemy w prosty sposób sprawdzić statystyki opisowe naszego zbioru danych.
summary(landsat)
## B1 B2 B3 B4
## Min. : 6577 Min. : 32 Min. : 6844 Min. : 7686
## 1st Qu.: 8039 1st Qu.: 8220 1st Qu.: 9093 1st Qu.: 8649
## Median : 8425 Median : 8614 Median : 9609 Median : 9353
## Mean : 8650 Mean : 8921 Mean : 9903 Mean : 9961
## 3rd Qu.: 9044 3rd Qu.: 9370 3rd Qu.:10396 3rd Qu.:10796
## Max. :19736 Max. :21314 Max. :24431 Max. :27674
## NA's :48391 NA's :48390 NA's :48390 NA's :48390
## B5 B6 B7
## Min. : 7240 Min. : 7378 Min. : 7357
## 1st Qu.:14669 1st Qu.:12482 1st Qu.:10002
## Median :17462 Median :14082 Median :11417
## Mean :18004 Mean :14673 Mean :12519
## 3rd Qu.:21254 3rd Qu.:16616 3rd Qu.:14170
## Max. :31081 Max. :50865 Max. :60130
## NA's :48390 NA's :48390 NA's :48390
Jak możemy zauważyć wartości odbicia spektralnego dla naszego zbioru danych są w zakresie od kilku do kilkunastu tysięcy dla każdego kanału. Odbicie spektralne powinno być w przedziale od 0 do 1, w związku z czym nasze dane musimy przeskalować za pomocą poniższego równania:
\[x = x \cdot 0.0000275 - 0.2\]
Dla przykładu, wartość piksela w kanale bliskiej podczerwieni wynosi 15000. Używając powyższego wzoru musimy tę wartość przemnożyć przez 0,0000275 (scale factor), a następnie odjąć 0,2 (offset). Jako wynik otrzymamy odbicie o wartości równej 0,2125. Należy pamiętać, że każdy produkt/kolekcja posiada inny wzór i konieczne jest zapoznanie się z dokumentacją.
Nie ma potrzeby stosowania tego wzoru osobno dla każdego kanału w pętli, ponieważ operacje matematyczne w pakiecie terra są domyślnie stosowane dla wszystkich kanałów.
landsat = landsat * 2.75e-05 - 0.2
summary(landsat)
## B1 B2 B3 B4
## Min. :-0.02 Min. :-0.20 Min. :-0.01 Min. :0.01
## 1st Qu.: 0.02 1st Qu.: 0.03 1st Qu.: 0.05 1st Qu.:0.04
## Median : 0.03 Median : 0.04 Median : 0.06 Median :0.06
## Mean : 0.04 Mean : 0.05 Mean : 0.07 Mean :0.07
## 3rd Qu.: 0.05 3rd Qu.: 0.06 3rd Qu.: 0.09 3rd Qu.:0.10
## Max. : 0.34 Max. : 0.39 Max. : 0.47 Max. :0.56
## NA's :48391 NA's :48390 NA's :48390 NA's :48390
## B5 B6 B7
## Min. :0.00 Min. :0.00 Min. :0.00
## 1st Qu.:0.20 1st Qu.:0.14 1st Qu.:0.08
## Median :0.28 Median :0.19 Median :0.11
## Mean :0.30 Mean :0.20 Mean :0.14
## 3rd Qu.:0.38 3rd Qu.:0.26 3rd Qu.:0.19
## Max. :0.65 Max. :1.20 Max. :1.45
## NA's :48390 NA's :48390 NA's :48390
Nadal możemy zauważyć, że pewne wartości przekraczają nasz zakres od 0 do 1. Są to wartości odstające, które zazwyczaj związane są z błędnym pomiarem lub nadmierną saturacją. Można ten problem rozwiązać na dwa sposoby:
NA).Pierwszy sposób może spowodować, że stracimy dużą część zbioru danych. Natomiast drugi sposób może powodować przekłamania.
# sposób nr 1
landsat[landsat < 0] = NA
landsat[landsat > 1] = NA
# sposób nr 2
landsat[landsat < 0] = 0
landsat[landsat > 1] = 1
Po przeskalowaniu wartości możemy wyświetlić kompozycję RGB. W tym
przypadku zamiast funkcji plot() należy użyć funkcji
plotRGB() oraz zdefiniować kolejność kanałów czerwonego,
zielonego oraz niebieskiego. Oprócz tego należy podać maksymalną wartość
odbicia dla kanałów (w naszym przypadku scale = 1). Często
zdarza się, że kompozycje są zbyt ciemne/jasne, wtedy warto zastosować
rozciągnięcie kolorów używając argumentu stretch = "lin"
lub stretch = "hist".
# plotRGB(landsat, r = 4, g = 3, b = 2, scale = 1)
plotRGB(landsat, r = 4, g = 3, b = 2, scale = 1, stretch = "lin")
# wczytanie pakietu do grupowania danych
library("cluster")
Dane do modelowania muszą zostać przygotowane w odpowiedni sposób.
Modele klasyfikacyjne najczęściej na etapie trenowania wymagają macierzy
lub ramki danych (data frame). Dane rastrowe można przetworzyć
do macierzy przy użyciu funkcji values().
mat = values(landsat)
nrow(mat) # wyświetla liczbę wierszy
## [1] 1238760
Za pomocą interaktywnej funkcji View() możemy sprawdzić
jak wygląda nasza macierz.
View(mat)
Jak można zauważyć, mnóstwo jej wartości oznaczonych jest jako brak
danych (głównie są to wartości poza obszarem analizy). Zazwyczaj modele
nie obsługują NA, więc musimy je usunąć. Służy do tego
dedykowana funkcja na.omit().
mat_omit = na.omit(mat)
nrow(mat_omit)
## [1] 640802
Teraz przejdziemy do najważniejszego etapu analizy, czyli do
wytrenowania modelu. Użyjemy prostego modelu grupowania metodą
k-średnich (k-means). Ten model wymaga jedynie, aby podać z
góry liczbę grup/klastrów (argument centers). Jest to
algorytm stochastyczny, więc za każdym razem zwraca inne wyniki. Żeby
analiza była powtarzalna musimy ustawić ziarno losowości –
set.seed().
set.seed(1)
mdl = kmeans(mat_omit, centers = 5)
W wyniku powyższej operacji otrzymaliśmy m.in.:
mdl$centers).mdl$cluster).Wyświetlmy te obiekty:
mdl$centers
## B1 B2 B3 B4 B5 B6 B7
## 1 0.02308493 0.02740529 0.04435163 0.04147431 0.1727330 0.1309803 0.07958614
## 2 0.02044635 0.02579641 0.05547638 0.03750665 0.4421821 0.1492823 0.07336897
## 3 0.08051378 0.09522015 0.13257100 0.16603768 0.2716112 0.3437177 0.32147881
## 4 0.03431384 0.04119050 0.07134604 0.06554308 0.3509450 0.2042120 0.12686408
## 5 0.04921918 0.05830302 0.08442430 0.09825275 0.2417769 0.2551757 0.19605182
head(mdl$cluster) # wyświetla pierwsze 6 elementów
## [1] 4 5 5 1 1 1
Oznacza to, że pierwszy wiersz (reprezentujący pojedyncze oczko siatki) należy do grupy 3, drugi do grupy 2, trzeci do grupy 2, itd. Kolejnym etapem jest stworzenie mapy na podstawie otrzymanego wektora z klastrami.
Na początku musimy przygotować pusty wektor składający się z
całkowitej liczby pikseli rastra. Można to sprawdzić za pomocą funkcji
ncell(). W naszym przypadku jest to 1 238 760.
vec = rep(NA, ncell(landsat)) # przygotuj pusty wektor
Następnie musimy przypisać nasze grupy w wektorze w odpowiednie
miejsca, tj. tym, które nie są zamaskowane (NA). Do
niezamaskowanych wartości można odwołać się przez funkcję
complete.cases().
# zastąp tylko te wartości, które nie są NA
vec[complete.cases(mat)] = mdl$cluster
W ostatnim kroku należy skopiować metadane obiektu
landsat, ale tylko z jedną warstwą, i przypisać mu wartości
wektora vec.
clustering = rast(landsat, nlyrs = 1, vals = vec)
Sprawdźmy teraz jak wyglądają utworzone grupy na mapie.
colors = rainbow(5, alpha = NULL) # wybierz 5 kolorów z wbudowanej palety `rainbow`
plot(clustering, type = "classes", col = colors)
Istotą grupowania jest stworzenie grupy składających się z podobnych elementów. Natomiast naszym zadaniem jest interpretacja, co przedstawiają utworzone grupy oraz ich nazwanie. Interpretacja jest trudnym zadaniem, a często jej wyniki są niejasne. W tym celu, niezbędna jest analiza statystyk opisowych grup oraz wspomaganie się różnymi kompozycjami (w naturalnych oraz fałszywszy kolorach). Bardzo przydatna jest również wiedza o właściwościach spektralnych obiektów.
Poniżej znajduje się przykładowy wynik takiej interpretacji.
colors = c("#086209", "#fdd327", "#d9d9d9", "#29a329", "#91632b")
category = c("lasy/woda", "pola uprawne", "odkryta gleba", "roślinność", "nieużytki")
plot(clustering, col = colors, type = "classes", levels = category)
Jeśli wynik jest satysfakcjonujący, to możemy go zapisać używając
funkcji writeRaster(). Taki plik można później wczytać w
R lub innym programie obsługującym dane przestrzenne
(np. QGIS).
writeRaster(clustering, "clustering.tif")
Do analizy i charakteryzacji grup, zamiast statystyk opisowych, mogą zostać wykorzystane wizualizacje. Największe możliwości dostarcza pakiet ggplot2. Tutaj można znaleźć darmowy podręcznik oraz gotowe “przepisy”.
ggplot2 wymaga przygotowania zbioru danych do odpowiedniej postaci. Dane muszą być przedstawione jako ramka danych w tzw. formie długiej (wiele wierszy), podczas gdy standardowe funkcje do wizualizacji wymagają formy szerokiej (wiele kolumn). Takiej zmiany można dokonać w prosty sposób używając pakietu tidyr.
# install.packages(c("tidyr", "ggplot2"))
library("tidyr") # transformacja danych
library("ggplot2") # wizualizacja danych
Nasz zbiór danych jest całkiem pokaźny (blisko 9 mln wartości),
jednak nie ma potrzeby przedstawiania wszystkich danych na wykresie.
Wymaga to więcej RAM i znacząco wydłuża czas rysowania. Prawie
identyczny efekt można uzyskać wykorzystując mniejszą próbkę danych.
Jako przykład zobrazujmy jedynie 10 000 wartości z każdego kanału
spektralnego. Do stworzenia losowej próby służy funkcja
sample(). W wyniku otrzymamy indeksy wylosowanych
wierszy.
idx = sample(1:nrow(mat_omit), size = 10000)
head(idx) # wyświetl 6 pierwszy indeksów
## [1] 294762 392686 640775 538191 270373 549593
Połączmy teraz wylosowane wiersze z macierzy z odpowiednimi grupami
(cbind()). Następnie macierz zamienimy na ramkę danych
(as.data.frame()).
stats = cbind(mat_omit[idx, ], cluster = mdl$cluster[idx])
stats = as.data.frame(stats)
head(stats)
## B1 B2 B3 B4 B5 B6 B7 cluster
## 1 0.0215950 0.0250875 0.0423575 0.0363350 0.1958350 0.1312650 0.0719200 1
## 2 0.0209625 0.0296525 0.0567400 0.0480775 0.3601475 0.1481500 0.0784650 4
## 3 0.0611400 0.0707100 0.1034350 0.1209800 0.3034425 0.3146350 0.2814150 3
## 4 0.0356200 0.0396350 0.0560250 0.0620475 0.1954225 0.2082925 0.1317050 1
## 5 0.0338875 0.0392775 0.0642200 0.0594900 0.3336650 0.2240500 0.1200450 4
## 6 0.0487650 0.0554750 0.0807475 0.0867425 0.2840550 0.2358200 0.1891525 5
Jak można zauważyć, powyższe dane mają formę szeroką (każdy kanał
spektralny zapisany jest w osobnej kolumnie). Teraz musimy zmienić
formę, w której otrzymamy dwie kolumn – kanał oraz wartość. W tym celu
wykorzystamy funkcję pivot_longer().
stats = pivot_longer(stats, cols = 1:7, names_to = "band", values_to = "value")
Dla formalności możemy jeszcze zmienić typ danych (klastrów i kanałów) na kategoryczny (factor). W praktyce związane jest to z uproszczeniem struktury danych (przejście ze skali ilorazowej do nominalnej).
stats$cluster = as.factor(stats$cluster)
stats$band = as.factor(stats$band)
head(stats)
## # A tibble: 6 × 3
## cluster band value
## <fct> <fct> <dbl>
## 1 1 B1 0.0216
## 2 1 B2 0.0251
## 3 1 B3 0.0424
## 4 1 B4 0.0363
## 5 1 B5 0.196
## 6 1 B6 0.131
Struktura danych jest już przygotowana. Teraz stwórzmy prosty wykres pudełkowy.
ggplot(stats, aes(x = band, y = value, fill = cluster)) +
geom_boxplot()
Zmieńmy kilka domyślnych parametrów żeby poprawić odbiór ryciny.
ggplot(stats, aes(x = band, y = value, fill = cluster)) +
geom_boxplot(show.legend = FALSE) +
scale_fill_manual(values = colors) +
facet_wrap(vars(cluster)) +
xlab("Kanał") +
ylab("Odbicie") +
theme_light()
Zmieniając facet_wrap(vars(cluster)) na
facet_wrap(vars(band)), zamiast zestawienia kanałów w
poszczególnych panelach, możemy zestawić grupy.